/****************************************************************************
 * Copyright (C) 2009-2010 SciTouch LLC
 * 
 * This file is part of Indigo toolkit.
 * 
 * This file may be distributed and/or modified under the terms of the
 * GNU General Public License version 3 as published by the Free Software
 * Foundation and appearing in the file LICENSE.GPL included in the
 * packaging of this file.
 * 
 * This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
 * WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
 ***************************************************************************/

#include "molecule/molecule_aromatic_cycles_extractor.h"

#include "molecule/molecule.h"
#include "graph/cycle_enumerator.h"

MoleculeAromaticCyclesExtractor::MoleculeAromaticCyclesExtractor (
   const Molecule &mol) 
   : _mol(mol),
   TL_CP_GET(_cycles_vertices),
   TL_CP_GET(_cycles_edges),
   TL_CP_GET(_cycles_offset),
   TL_CP_GET(_aromatic_groups),
   TL_CP_GET(_bond_aromatic_group)
{
   _cycles_vertices.clear();
   _aromatic_groups.clear();
   _bond_aromatic_group.clear_resize(_mol.edgeEnd());
   _bond_aromatic_group.fffill();
}

void MoleculeAromaticCyclesExtractor::extractBondAromaticCycles (
   int bond, Array<int> &cycles_indices, int max_cycle_length)
{
   cycles_indices.clear();

   int bond_type  = _mol.getBond(bond).type;
   if (bond_type != BOND_AROMATIC)
      return;
   
   if (_bond_aromatic_group[bond] == -1)
   {
      // Cycles were not collected
      _collectCycles(bond, max_cycle_length);
   }

   int group_idx = _bond_aromatic_group[bond];
   const AromaticGroup &group = _aromatic_groups[group_idx];
   if (group.max_cycle_length < max_cycle_length)
      throw Error("max_cycle_length should be the same for "
                  "every request for specific aromatic system");

   for (int i = group.start_cycles; i < group.end_cycles; i++)
   {
      if (_hasCycleBond(i, bond))
         cycles_indices.push(i);
   }
}

bool MoleculeAromaticCyclesExtractor::_hasCycleBond  (int cycle_idx, int bond) const
{
   int len = getCycleLength(cycle_idx);
   const int *edges = getCycleEdges(cycle_idx);
   for (int i = 0; i < len; i++)
      if (edges[i] == bond)
         return true;

   return false;
}

int MoleculeAromaticCyclesExtractor::getCycleLength (int cycle_idx) const
{
   int offset = _cycles_offset[cycle_idx];
   if (cycle_idx < _cycles_offset.size() - 1)
      return _cycles_offset[cycle_idx + 1] - offset;

   return _cycles_vertices.size() - offset;
}

const int* MoleculeAromaticCyclesExtractor::getCycleVertices (int idx) const
{
   int offset = _cycles_offset[idx];
   return _cycles_vertices.ptr() + offset;
}

const int* MoleculeAromaticCyclesExtractor::getCycleEdges (int idx) const
{
   int offset = _cycles_offset[idx];
   return _cycles_edges.ptr() + offset;
}

void MoleculeAromaticCyclesExtractor::_extractAromaticSystem (int bond, 
      Array<int> &vertices, Array<int> &edges) const
{
   // Find aromatic system
   enum { STATE_FREE = 0, STATE_IN_QUEUE, STATE_PROCESSED };
   QS_DEF(Array<int>, vertices_state);
   vertices_state.clear_resize(_mol.vertexEnd());
   vertices_state.zerofill();

   QS_DEF(Array<int>, queue);
   const Edge &edge = _mol.getEdge(bond);
   queue.push(edge.beg);

   vertices.clear();
   edges.clear();

   vertices.push(edge.beg);
   vertices_state[edge.beg] = STATE_IN_QUEUE;

   while (queue.size() != 0)
   {
      int v = queue.top();
      queue.pop();

      const Vertex &vertex = _mol.getVertex(v);
      for (int nei = vertex.neiBegin(); 
               nei != vertex.neiEnd(); 
               nei = vertex.neiNext(nei))
      {
         int nei_vertex = vertex.neiVertex(nei);
         if (vertices_state[nei_vertex] == STATE_PROCESSED)
            continue;

         int nei_edge = vertex.neiEdge(nei);
         int nei_bond_type = _mol.getBond(nei_edge).type;
         if (nei_bond_type != BOND_AROMATIC)
            continue;

         edges.push(nei_edge);

         if (vertices_state[nei_vertex] != STATE_IN_QUEUE)
         {
            queue.push(nei_vertex);

            vertices.push(nei_vertex);

            vertices_state[nei_vertex] = STATE_IN_QUEUE;
         }
      }
      vertices_state[v] = STATE_PROCESSED;
   }
}

void MoleculeAromaticCyclesExtractor::_extractPiLables (const Array<int> &vertices, 
         Array<int> &pi_labels) const
{
   pi_labels.clear_resize(vertices.size());

   for (int i = 0; i < vertices.size(); i++)
   {
      int v_idx = vertices[i];

      int min_connectivity = 0;
      const Vertex &vertex = _mol.getVertex(v_idx);
      for (int nei = vertex.neiBegin(); 
               nei != vertex.neiEnd(); 
               nei = vertex.neiNext(nei))
      {
         int nei_edge_idx = vertex.neiEdge(nei);
         int bond_type = _mol.getBond(nei_edge_idx).type;

         if (bond_type == BOND_AROMATIC)
            min_connectivity++;
         else 
            min_connectivity += bond_type;
      }

      int valence, hyd;
      const Atom &atom = _mol.getAtom(v_idx);
      Element::calcValence(atom.label, atom.charge, atom.radical, 
         min_connectivity, valence, hyd, true);

      if (hyd != atom.implicit_h)
      {
         // This atom was drawn with double bound
         pi_labels[i] = 1;
      }
      else
      {
         // Double bond wasn't attached
         int lonepairs;
         int vac = _mol.getVacantPiOrbitals(v_idx, min_connectivity, &lonepairs);
         if (vac > 0)
            pi_labels[i] = 0;
         else if (lonepairs > 0)
            pi_labels[i] = 2;
         else
            throw Error("Unable to determine atom pi-label");
      }
   }
}

bool MoleculeAromaticCyclesExtractor::_handleCycle (
                  const Graph &graph, const Array<int> &vertices, 
                  const Array<int> &edges, void *context_ptr)
{
   EnumerationContext &context = *(EnumerationContext *)context_ptr;
   MoleculeAromaticCyclesExtractor &self = *context.self;

   // Check pi-lables sum
   int sum = 0;
   for (int i = 0; i < vertices.size(); i++)
   {
      int vertex_in_system = vertices[i];
      sum += context.pi_labels[vertex_in_system];
   }

   if (sum % 4 != 2)
      return true;

   if (vertices.size() != edges.size())
      throw Error("Internal logic error");

   // Cycle is aromatic
   self._cycles_offset.push(self._cycles_vertices.size());
   for (int i = 0; i < vertices.size(); i++)
   {
      int v = context.aromatic_vertices[vertices[i]];
      int e = context.aromatic_edges[edges[i]];
      self._cycles_vertices.push(v);
      self._cycles_edges.push(e);
   }
   return true;
}


void MoleculeAromaticCyclesExtractor::_collectCycles (int bond, int max_cycle_length)
{
   QS_DEF(EnumerationContext, context);
   context.self = this;

   _extractAromaticSystem(bond, context.aromatic_vertices, context.aromatic_edges);

   _extractPiLables(context.aromatic_vertices, context.pi_labels);

   QS_DEF(Graph, aromatic_system);
   aromatic_system.makeEdgeSubgraph(_mol, context.aromatic_vertices, 
      context.aromatic_edges, NULL, NULL);

   AromaticGroup &group = _aromatic_groups.push();
   group.start_cycles = _cycles_offset.size();
   group.max_cycle_length = max_cycle_length;

   // Enumrate cycles
   CycleEnumerator enumerator(aromatic_system);
   if (max_cycle_length != -1)
      enumerator.max_length = max_cycle_length;
   enumerator.context = &context;
   enumerator.cb_handle_cycle = _handleCycle;
   enumerator.process();

   group.end_cycles = _cycles_offset.size();

   // Mark bond
   for (int i = 0; i < context.aromatic_edges.size(); i++)
   {
      int edge = context.aromatic_edges[i];
      _bond_aromatic_group[edge] = _aromatic_groups.size() - 1;
   }
}
